*! version 1.0.0 13July 2017

program randcmd, eclass
	version 13.1
	syntax anything , treatvars(string) [calc1(string) calc2(string) calc3(string) calc4(string) calc5(string) calc6(string) calc7(string) calc8(string) calc9(string) calc10(string) reps(integer 1000) strata(string) groupvar(string) seed(integer 1) saving(string)]
	tempname b f bb ff fff info T ResB ResSE ResF ResCoef ResEqn ResOmni cov 
	tempvar U Order M n

preserve

	local oldseed = "`c(seed)'"

*Extracting list of treatment variables and post-treatment calculations, establishing sample (treatvars ~= .), checking treatvars do not vary within groupvar
	unab treatvars: `treatvars'
	local calc = 0
	forvalues k = 1/10 {
		if ("`calc`k''" ~= "") local calc = `k'
		}
	foreach var in `treatvars' {
		quietly drop if `var' == .
		}
	local error = 0
	if ("`groupvar'" ~= "") {
		foreach var in `treatvars' {
			quietly egen `M' = sd(`var'), by(`groupvar')
			quietly sum `M'
			if (r(mean) ~= 0) {
				display as error "`var' varies within `groupvar'.  Base treatment variables should not vary within treatment groupings."
				local error = 1
				}
			quietly drop `M'
			}
		}
	if (`error' == 1) exit
			

*Baseline estimating equations
	local eqn = 0
	local treatnumber = 0
	while "`anything'" ~= "" {
		local eqn = `eqn' + 1
		gettoken eqn`eqn' anything: anything, match(match)
		gettoken treat`eqn' eqn`eqn': eqn`eqn', match(match)
		local treatnumber = `treatnumber' + wordcount("`treat`eqn''")
		if (wordcount("`treat`eqn''") == 0) {
			display as error "No treatment variables specified for equation `eqn'."
			local error = 1
			}		
		}
	if (`error' == 1) exit	
	matrix `f' = J(`eqn',3,.)
	matrix `b' = J(`treatnumber',2,.)
	local c = 0
	local length = 0
	forvalues k = 1/`eqn' {
		`eqn`k''
		local test`k'  = ""
		foreach var in `treat`k'' {
			local length = max(`length',length("`var'"))
			local c = `c' + 1
			capture matrix `b'[`c',1] = _b[`var'], _se[`var']
			if (_rc ~= 0) {
				display as error "`var' cannot be found in estimating equation `k'."
				local error = 1
				}
			local test`k' = "`test`k''" + "(_b[`var']==0)"
			}
		if (`error' == 1) continue, break
		test `test`k''
		local f`k' = r(df)
		if (r(df_r) ~= .) local ftype`k' = "r(F)"
		if (r(df_r) == .) local ftype`k' = "r(chi2)"
		matrix `f'[`k',1] = `ftype`k'', r(p), `f`k''
		display " "
		}
	if (`error' == 1) exit
	matrix `info' = J(`eqn',2,1)
	forvalues k = 1/`eqn' {
		matrix `info'[`k',2] = `info'[`k',1] + wordcount("`treat`k''") - 1
		if (`k' < `eqn') matrix `info'[`k'+1,1] = `info'[`k',2] + 1
		}


*Displaying treatment variables so that user can confirm that programme has correctly identified treatment variables and interaction equations
	display "Treatment variables determined directly by randomization: `treatvars'.", _newline
	display "Post-randomization treatment based calculations:  `calc'."
	forvalues k = 1/`calc' {
		display "  `k':   `calc`k''"
		}
	display " "
	forvalues k = 1/`eqn' {
		display "Treatment based variables tested in equation `k': `treat`k''."
		if  ("`treat`k''" == "") {
			display as error "No treatment variables in equation `k'."
			local error = 1
			}
		}
	if (`error' == 1) exit
	display " "


*Testing Calculations
	if (`calc' > 0) {
		display "Testing calculations.  Should execute without error and, given that base treatment has not been rerandomized, produce unchanged treatment coefficients."
		forvalues k = 1/`calc' {
			display "  `k':   `calc`k''"
			quietly `calc`k''
			}
		matrix `bb' = J(`treatnumber',2,.)
		forvalues k = 1/`eqn' {
			capture `eqn`k''
			if (_rc == 0) {
				local c = `info'[`k',1]
				foreach var in `treat`k'' {
					capture matrix `bb'[`c',1] = _b[`var'], _se[`var']
					local c = `c' + 1
					}
				}
			matrix `ff' = `b'[`info'[`k',1]..`info'[`k',2],1..2] - `bb'[`info'[`k',1]..`info'[`k',2],1..2]
			matrix `ff' = trace(`ff''*`ff')
			if (`ff'[1,1] > 1e-06) {
				display as error "Coefficients in equation `k' have changed.  Results of calculations are inconsistent with current data."
				local error = 1
				}
			}
		}
	if (`error' == 1) exit


*Preparing variables and matrices to be used in randomization analysis
	set seed `seed'
	if ("`groupvar'" ~= "") {
		egen `M' = group(`groupvar')
		quietly sum `M'
		local N = r(max)
		quietly bysort `M': gen `n' = _n
		sort `n' `strata' `M'
		quietly generate `Order' = _n
		}
	if ("`groupvar'" == "") {
		local N = _N
		quietly generate `Order' = _n	
		sort `strata' `Order'
		}
	quietly generate double `U' = .
	mata list = J(1,0,"")
	foreach var in `treatvars' {
		mata list = list, "`var'"
		}
	mata `T' = st_data((1,`N'),list); `ResB' = J(`reps',`treatnumber',.); `ResSE' = J(`reps',`treatnumber',.); `ResF' = J(`reps',`eqn',.)

display " "
display "Running `reps' randomization iterations:"
		

*Randomization iterations
	forvalues count = 1/`reps' {
		if (ceil(`count'/10)*10 == `count') display "`count'", _continue

*Randomizing direct treatment and recalculating treatment based variables
		if ("`groupvar'" == "") {
			quietly sort `strata' `Order'
			quietly replace `U' = uniform()
			quietly sort `strata' `U'	
			mata st_store(.,list,`T')
			}
		if ("`groupvar'" ~= "") {
			quietly sort `n' `strata' `Order'  
			quietly replace `U' = uniform() if _n <= `N'
			quietly sort `strata' `U' in 1/`N'
			mata st_store((1,`N'),list,`T')
			quietly sort `M' `n'
			foreach var in `treatvars' {
				quietly replace `var' = `var'[_n-1] if `n' > 1
				}
			}						
		forvalues k = 1/`calc' {
			quietly `calc`k''
			}

*Estimating equations
		matrix `bb' = J(`treatnumber',2,.)
		matrix `ff' = J(`eqn',1,.)
		forvalues k = 1/`eqn' {
			capture `eqn`k''
			if (_rc == 0) {
				local c = `info'[`k',1]
				foreach var in `treat`k'' {
					capture matrix `bb'[`c',1] = _b[`var'], _se[`var']
					local c = `c' + 1
					}
				capture test `test`k''
				if (_rc == 0 & r(drop) == 0 & r(df) == `f`k'') capture matrix `ff'[`k',1] = `ftype`k''
				}
			}
		mata `bb' = st_matrix("`bb'"); `ff' = st_matrix("`ff'"); `ResB'[`count',1...] = `bb'[1...,1]'; `ResSE'[`count',1...] = `bb'[1...,2]'; `ResF'[`count',1...] = `ff'[1...,1]'
		}

display, _newline


*Calculating p-values
	mata `b' = st_matrix("`b'"); `f' = st_matrix("`f'"); `ResCoef' = J(`treatnumber',6,.); `ResEqn' = J(`eqn',6,.)
	forvalues c = 1/`treatnumber' {
		mata `bb' = (`ResB'[1...,`c']:~=.); `bb' = `bb':*(`ResSE'[1...,`c']:~=.); `bb' = `bb':*(`ResSE'[1...,`c']:~=0)
		mata `ff' = select(`ResB'[1...,`c'],`bb'); `ff' = (abs(`ff'):>abs(`b'[`c',1])*1.000001), (abs(`ff'):>abs(`b'[`c',1])*.999999); `ff' = colsum(`ff'), rows(`ff')
		mata `ResCoef'[`c',1..3] = `ff'[1,1]/(`ff'[1,3]+1), (`ff'[1,2]+1)/(`ff'[1,3]+1), `ff'[1,3]
		mata `ff' = select(`ResB'[1...,`c']:/`ResSE'[1...,`c'],`bb'); `ff' = (abs(`ff'):>abs(`b'[`c',1]/`b'[`c',2])*1.000001), (abs(`ff'):>abs(`b'[`c',1]/`b'[`c',2])*.999999); `ff' = colsum(`ff'), rows(`ff')
		mata `ResCoef'[`c',4..6] = `ff'[1,1]/(`ff'[1,3]+1), (`ff'[1,2]+1)/(`ff'[1,3]+1), `ff'[1,3]
		}
	forvalues e = 1/`eqn' {
		mata `bb' = (`ResF'[1...,`e']:~=.)
		local a1 = `info'[`e',1]
		local a2 = `info'[`e',2]
		forvalues c = `a1'/`a2' {
			mata `bb' = `bb':*(`ResB'[1...,`c']:~=.); `bb' = `bb':*(`ResSE'[1...,`c']:~=.); `bb' = `bb':*(`ResSE'[1...,`c']:~=0)
			}	
		mata `ff' = select(`ResB'[1...,`a1'..`a2'],`bb'); `cov' = `ff':-mean(`ff'); `cov' = `cov''*`cov'/rows(`cov'); `cov' = invsym(`cov')
		mata `ff' = rowsum(`ff'*`cov':*`ff'); `fff' = `b'[`a1'..`a2',1]'*`cov'*`b'[`a1'..`a2',1]
		mata `ff' = (`ff':>`fff'*1.000001), (`ff':>`fff'*.999999); `ff' = colsum(`ff'), rows(`ff')
		mata `ResEqn'[`e',1..3] = `ff'[1,1]/(`ff'[1,3]+1), (`ff'[1,2]+1)/(`ff'[1,3]+1), `ff'[1,3]
		mata `ff' = select(`ResF'[1...,`e'],`bb'); `ff' = (`ff':>`f'[`e',1]*1.000001), (`ff':>`f'[`e',1]*.999999); `ff' = colsum(`ff'), rows(`ff')
		mata `ResEqn'[`e',4..6] = `ff'[1,1]/(`ff'[1,3]+1), (`ff'[1,2]+1)/(`ff'[1,3]+1), `ff'[1,3]
		}
	mata `bb' = J(`reps',1,1)
	forvalues c = 1/`treatnumber' {
		mata `bb' = `bb':*(`ResB'[1...,`c']:~=.); `bb' = `bb':*(`ResSE'[1...,`c']:~=.); `bb' = `bb':*(`ResSE'[1...,`c']:~=0)
		}
	mata `ff' = select(`ResB',`bb'); `cov' = `ff':-mean(`ff'); `cov' = `cov''*`cov'/rows(`cov'); `cov' = invsym(`cov')
	mata `ff' = rowsum(`ff'*`cov':*`ff'); `fff' = `b'[1...,1]'*`cov'*`b'[1...,1]
	mata `ff' = (`ff':>`fff'*1.000001), (`ff':>`fff'*.999999); `ff' = colsum(`ff'), rows(`ff')
	mata `ResOmni' = `ff'[1,1]/(`ff'[1,3]+1), (`ff'[1,2]+1)/(`ff'[1,3]+1), `ff'[1,3]
	mata `bb' = uniform(rows(`ResCoef'),1); `ResCoef' = `ResCoef'[1...,1..2], `ResCoef'[1...,1]+`bb':*(`ResCoef'[1...,2]-`ResCoef'[1...,1]), `ResCoef'[1...,4..5], `ResCoef'[1...,4]+`bb':*(`ResCoef'[1...,5]-`ResCoef'[1...,4]), `ResCoef'[1...,6]
	mata `bb' = uniform(rows(`ResEqn'),1); `ResEqn' = `ResEqn'[1...,1..2], `ResEqn'[1...,1]+`bb':*(`ResEqn'[1...,2]-`ResEqn'[1...,1]), `ResEqn'[1...,4..5], `ResEqn'[1...,4]+`bb':*(`ResEqn'[1...,5]-`ResEqn'[1...,4]), `ResEqn'[1...,6]
	mata `bb' = uniform(rows(`ResOmni'),1); `ResOmni' = `ResOmni'[1...,1..2], `ResOmni'[1...,1]+`bb':*(`ResOmni'[1...,2]-`ResOmni'[1...,1]), `ResOmni'[1...,3]
	mata st_matrix("`ResCoef'",`ResCoef'); st_matrix("`ResEqn'",`ResEqn'); st_matrix("`ResOmni'",`ResOmni')


*Displaying results
local length = `length' + 3 + (`eqn' > 9)
local length = max(`length',9)
local a1 = `length' + 6
forvalues k = 2/7 {
	local a`k' = `a1' + 13*(`k'-1)
	}
local aa1 = `a1' - 4
local aa2 = floor(`length'/2)-1
local aa3 = `a2' - 3
local aa5 = `a5' - 3

display as text " "
display "Randomization p-values for individual coefficients:", _newline
display as text _col(`aa3') %15s  "randomization-c" _col(`aa5') "randomization-t" 
display as text _col(`aa2') %8s "equation:" _col(`a1') %8s "minimum" _col(`a2') %8s "maximum" _col(`a3') %10s "randomized" _col(`a4') %8s "minimum" _col(`a5') %8s "maximum" _col(`a6') %10s "randomized" _col(`a7') %10s "successful"
display as text _col(`aa2') %8s "variable" _col(`a1') %8s "p-value" _col(`a2') %8s "p_value" _col(`a3') %8s "p-value" _col(`a4') %8s "p-value" _col(`a5') %8s "p-value" _col(`a6') %8s "p-value" _col(`a7') %10s "iterations"
display as text "{hline `aa1'}{c +}{hline 90}"
	local i = 0
	forvalues e = 1/`eqn' {
		foreach var in `treat`e'' {
			local i = `i' + 1
			display as text _col(2) %-`length's "`e': `var'" _col(`aa1') " {c |}" , _continue
			display as result _col(`a1') %7.6g `ResCoef'[`i',1] _col(`a2') %7.6g `ResCoef'[`i',2] _col(`a3') %7.6g `ResCoef'[`i',3] _col(`a4') %7.6g `ResCoef'[`i',4] _col(`a5') %7.6g `ResCoef'[`i',5] _col(`a6') %7.6g `ResCoef'[`i',6] _col(`a7') %7.6g `ResCoef'[`i',7]  
			}
		}
display as text "{hline `aa1'}{c BT}{hline 90}", _newline

display " "
display "Randomization p-values for equation-level joint test of treatment significance:", _newline
display as text _col(24) %15s  "randomization-c" _col(63) "randomization-t" 
display as text _col(14) %8s "minimum" _col(27) %8s "maximum" _col(39) %10s "randomized" _col(53) %8s "minimum" _col(66) %8s "maximum" _col(78) %7s "randomized" _col(93) %10s "successful"
display as text _col(2) %8s "equation" _col(14) %8s "p-value" _col(27) %8s "p_value" _col(41) %7s "p-value" _col(53) %8s "p-value" _col(66) %8s "p-value" _col(80) %7s "p-value" _col(93) %10s "iterations"
display as text "{hline 9}{c +}{hline 92}"
	forvalues e = 1/`eqn' {
		display as text _col(5) %-12s `e' _col(10) "{c |}" , _continue
		display as result _col(14) %7.6g `ResEqn'[`e',1] _col(27) %7.6g `ResEqn'[`e',2] _col(40) %7.6g `ResEqn'[`e',3] _col(53) %7.6g `ResEqn'[`e',4] _col(66) %7.6g `ResEqn'[`e',5] _col(79) %7.6g `ResEqn'[`e',6] _col(93) %7.6g `ResEqn'[`e',7]  
		}
display as text "{hline 9}{c BT}{hline 92}", _newline

display " "
display "Randomization p-value for omnibus joint test of treatment significance:", _newline
display as text _col(18) %15s  "randomization-c" 
display as text _col(2) %8s "minimum" _col(15) %8s "maximum" _col(27) %10s "randomized" _col(41) %10s "successful"
display as text _col(2) %8s "p-value" _col(15) %8s "p_value" _col(29) %7s "p-value" _col(41) %10s "iterations"
display as text _col(3) "{hline 48}"
	display as result _col(2) %7.6g `ResOmni'[1,1] _col(15) %7.6g `ResOmni'[1,2] _col(28) %7.6g `ResOmni'[1,3] _col(41) %7.6g `ResOmni'[1,4]
display as text _col(3) "{hline 48}", _newline

*ereturn matrices
matrix colnames `ResCoef' = "min-c pvalue" "max-c pvalue" "rand-c pvalue" "min-t pvalue" "max-t pvalue" "rand-t pvalue" "iterations"
matrix colnames `ResEqn' = "min-c pvalue" "max-c pvalue" "rand-c pvalue" "min-t pvalue" "max-t pvalue" "rand-t pvalue" "iterations"
matrix colnames `ResOmni' = "min-c pvalue" "max-c pvalue" "rand-c pvalue" "iterations"
local list = ""
forvalues e = 1/`eqn' {
	foreach var in `treat`e'' {
		local list = "`list'" + " `e':`var' "
		}
	}
matrix rownames `ResCoef' = `list'
local list = ""
forvalues e = 1/`eqn' {
	local list = "`list'" + " equation:`e' "
	}
matrix rownames `ResEqn' = `list'
matrix rownames `ResOmni' = "Omnibus test"

ereturn matrix RCoef = `ResCoef', copy
ereturn matrix REqn = `ResEqn', copy
ereturn matrix ROmni = `ResOmni', copy

set seed `oldseed'

*Saving, if user requested
if ("`saving'" ~= "") {
	quietly drop _all
	quietly set obs `reps'
	forvalues k = 1/`treatnumber' {
		quietly generate double ResB`k' = .
		}
	forvalues k = 1/`treatnumber' {
		quietly generate double ResSE`k' = .
		}
	forvalues k = 1/`eqn' {
		quietly generate double ResF`k' = .
		}
	mata st_store(.,.,(`ResB', `ResSE', `ResF'))
	save `saving'
	}

restore

end
